********************************************************************************
*********** This program is written for the present value project **************
********************************************************************************
* 05. Cross-sectional MB Decomposition using Portfolios May2021
********************************************************************************


********************************************************************************
* 05.00 Initialize STATA
********************************************************************************
clear all
set more off, perm

********************************************************************************
* 05.01  Settings
******************************************************************************** 
* Parameters for duration Calculation following Weber 2018
local discount_rate = 0.12 // discount rate 
local AR_ROE = 0.4067106129 // AR(1) coefficient of ROE, for terminal value
local ROE_SS = 0.12 // steady-state ROE, for terminal value
local AR_SG = 0.2411083953 // AR(1) coefficient BE equity growth
local SG_SS = 0.06 // steady-state sales growth rate
local horizon = 15

* Include Intangible
local intangible = 0

* Sample period
local start = 1965
local end = 2022

* Choose rho
local rho = 0.96 


********************************************************************************
* 05.02  Apply the sample period
********************************************************************************
* Load annual dataset
use "Data/Created/Annual Data 2021 with Bins", clear 

tempfile data
save "`data'"

keep bin_ME_L1 bin_bm_L1 year
drop if [_n]>=1
tempfile tmp
save "`tmp'"

********************************************************************************
* 05.03  Generate portfolio-level variables (time 0)
********************************************************************************  
use "`data'", clear

* Flag observations with missing firm-level variables
qui gen miss = 1 if mi(BE_L1) | mi(ME_L1) | mi(R) | mi(D) | mi(N_dilute) | ///
	mi(ME) | mi(BE) | mi(Y)

* Generate variables to be summed across firms within a portfolio
qui gen ME_L1_times_R_0 = ME_L1_times_R if miss!=1
qui gen ME_L1_0 = ME_L1 if miss!=1
qui gen BE_L1_0 = BE_L1 if miss!=1
qui gen ME_times_N_dilute_0 = ME_times_N_dilute if miss!=1
qui gen BE_L1_Y_0 = BE_L1 + Y if miss!=1
qui gen D_0 = D if miss!=1
qui gen BE_0 = BE if miss!=1
qui gen ME_0 = ME if miss!=1

* Current variables other than MB  
bys firm (year): g BE_L2 = BE_L1[_n-1] if year==year[_n-1]+1 & miss!=1
g BE_N_dilute = N_dilute * BE if year==year[_n-1]+1 & miss!=1

foreach var of varlist Y D BE_N_dilute N_dilute{
	bys firm (year): g `var'_L1 = `var'[_n-1] if year==year[_n-1]+1 & miss!=1 & ~mi(BE_L2)
	}

bys firm (year): g at_L1 = at[_n-1] if year==year[_n-1]+1 & miss!= 1 
bys firm (year): g at_L2 = at[_n-2] if year==year[_n-2]+2 & ~mi(at_L1) & miss!= 1 	 

g R_L1_num = R_L1 * ME_L1 if miss!= 1 
g R_L1_den = ME_L1 if ~mi(R_L1) & miss!= 1 

g ME_times_N_dilute_L1_0 = ME_L1*N_dilute_L1 if miss!=1

* Sum across firms within a portfolio
collapse (sum) ME_* BE_* D_* Y_L1 N_dilute_L1 at sale sale_L1 at_* R_L1_* ///
	, by(year binq_ME_L1 binq_bm_L1)

* Compute portfolio-level variables	
qui gen MB_L1_0 = ME_L1_0 / BE_L1_0
qui gen R_0 = ME_L1_times_R_0 / ME_L1_0
qui gen N_dilute_0 = ME_times_N_dilute_0 / ME_0
qui gen N_dilute_L1_0 = ME_times_N_dilute_L1_0 / ME_L1_0
qui gen ROE_0 = BE_L1_Y_0 / BE_L1_0 - 1
qui gen IVA_0 = (D_0 + N_dilute_0 * BE_0) / BE_L1_Y_0 - 1
qui gen MB_0 = ME_0 / BE_0
qui gen PVA_0 = 1 + (MB_0 - 1) * (N_dilute_0 * BE_0) / (D_0 + N_dilute_0 * BE_0)
qui gen sale_g = log(sale/sale_L1)	

* Drop if bins missing
drop if mi(binq_ME_L1) | mi(binq_bm_L1)
	
* Check that the exact nonlinear identity holds
qui gen R_implied_0 = (1+ROE_0)*(1+IVA_0)*PVA_0/MB_L1_0 - 1
qui gen err1_0 = R_0 - R_implied_0
summarize err1_0, detail	

* Generate log variables
gen mb_L1_0 = log(MB_L1_0)
gen r_0 = log(1+R_0)
gen roe_0 = log(1+ROE_0)
gen iva_0 = log(1+IVA_0)
gen mb_0 = log(MB_0)
gen pva_0 = log(PVA_0) - `rho'*mb_0 
gen ag_0 = log(at/at_L1)

* Current values
g r_L1_0 = log(1 + R_L1_num/R_L1_den)
g roe_L1_0 = log(1 + Y_L1/BE_L2)
g iva_L1_0 = log((D_L1 + BE_N_dilute_L1) / (BE_L2 + Y_L1))
g ag_L1_0 = log(at_L1/at_L2)	
g sgrowth_L1_0 = log(1/N_dilute_L1_0)   					// previously this variable was mislabeled as sgrowth_0
gen GHpf = 0
replace GHpf = 1 if sgrowth_L1_0 > 0.1					// portfolio analog of GH variable
replace GHpf = -1 if sgrowth_L1_0 < -0.005	
																
* Check the approximate linear identity (without pva)
qui gen r_implied_0 = roe_0 + iva_0 + `rho'*mb_0 - mb_L1_0  
qui gen err2_0 = r_0 - r_implied_0
summarize err2_0, detail	

* Append to temporary file and save
append using "`tmp'"
tempfile tmp
save "`tmp'"


********************************************************************************
* 05.04  Generate portfolio-level variables (time 1-19)
********************************************************************************
forvalues i = 0/19 {  

	use "`data'", clear
	local j = `i' - 1
	
	* Flag observations with missing firm-level variables 
	qui bys firm (year): gen miss = 1 if mi(BE[_n+`j']) | mi(ME[_n+`j']) | ///
		mi(R[_n+`i']) | mi(D[_n+`i']) | mi(N_dilute[_n+`i']) | ///
		mi(ME[_n+`i']) | mi(BE[_n+`i']) | mi(Y[_n+`i']) | year!=year[_n+`i']-`i'

	* Generate variables to be summed across firms within a portfolio
	qui bys firm (year): gen ME_L1_times_R_`i' = ME[_n+`j'] * R[_n+`i'] if miss!=1
	qui bys firm (year): gen ME_L1_`i' = ME[_n+`j'] if miss!=1
	qui bys firm (year): gen BE_L1_`i' = BE[_n+`j'] if miss!=1
	qui bys firm (year): gen ME_times_N_dilute_`i' = ME_times_N_dilute[_n+`i'] if miss!=1
	qui bys firm (year): gen BE_L1_Y_`i' = BE[_n+`j'] + Y[_n+`i'] if miss!=1
	qui bys firm (year): gen D_`i' = D[_n+`i'] if miss!=1
	qui bys firm (year): gen BE_`i' = BE[_n+`i'] if miss!=1
	qui bys firm (year): gen ME_`i' = ME[_n+`i'] if miss!=1
	qui bys firm (year): gen at_`i' = at[_n+`i'] if miss!=1
	qui bys firm (year): gen at_L1_`i' = at[_n+`j'] if miss!=1

	* Sum across firms within a portfolio
	collapse (sum) ME_* BE_* D_* at_*, by(year binq_ME_L1 binq_bm_L1)
		
	* Compute portfolio-level variables	
	qui gen R_`i' = ME_L1_times_R_`i' / ME_L1_`i'
	qui gen MB_L1_`i' = ME_L1_`i' / BE_L1_`i'
	qui gen N_dilute_`i' = ME_times_N_dilute_`i' / ME_`i'
	qui gen ROE_`i' = BE_L1_Y_`i' / BE_L1_`i' - 1
	qui gen IVA_`i' = (D_`i' + N_dilute_`i' * BE_`i') / BE_L1_Y_`i' - 1
	qui gen MB_`i' = ME_`i' / BE_`i'
	qui gen PVA_`i' = 1 + (MB_`i' - 1) * (N_dilute_`i' * BE_`i') / (D_`i' + N_dilute_`i' * BE_`i')	 

	* Drop if bins missing
	drop if mi(binq_ME_L1) | mi(binq_bm_L1)

	* Check that the exact nonlinear identity holds		
	gen R_implied_`i' = (1+ROE_`i')*(1+IVA_`i')*PVA_`i'/MB_L1_`i' - 1
	gen err1_`i' = R_`i' - R_implied_`i'
	summarize err1_`i', detail		

	* Generate log variables
	gen mb_L1_`i' = log(MB_L1_`i')
	gen r_`i' = log(1+R_`i')
	gen roe_`i' = log(1+ROE_`i')
	gen iva_`i' = log(1+IVA_`i')
	gen mb_`i' = log(MB_`i')
	gen pva_`i' = log(PVA_`i') - `rho'*mb_`i'
	gen ag_`i' = log(at_`i'/at_L1_`i')

	* Check the approximate linear identity (without pva)
	qui gen r_implied_`i' = roe_`i' + iva_`i' + `rho'*mb_`i' - mb_L1_`i'  
	qui gen err2_`i' = r_`i' - r_implied_`i'
	*summarize err2_`i', detail	
	
	* Merge with temporary file and save
	merge 1:1 year binq_ME_L1 binq_bm_L1 using "`tmp'", nogen
	tempfile tmp
	save "`tmp'"	
} 

* Compute portfolio-level duration, defined as in Equation (2) of Weber (2018, JFE): two versions, realized (consistent with CPV) and projected (as in Weber)
gen s_1 = (1 - `AR_SG') * `SG_SS' + `AR_SG' * sale_g
gen RonE_1 = (1 - `AR_ROE') * `ROE_SS' + `AR_ROE' * ROE_0
gen BV_1 = BE_0 * (1 + s_1)
gen CF_MW_1 = BE_0 * (1 + RonE_1) - BV_1
gen PV_CF_MW = CF_MW_1 / (1 + `discount_rate')
gen t_PV_CF_MW = CF_MW_1 / (1 + `discount_rate')

gen s_L1 = s_1
gen RonE_L1 = RonE_1
gen BV_L1 = BV_1

forvalues i = 2/`horizon' {
	gen s_`i' = (1 - `AR_SG') * `SG_SS' + `AR_SG' * s_L1
	gen RonE_`i' = (1 - `AR_ROE') * `ROE_SS' + `AR_ROE' * RonE_L1 
	gen BV_`i' = BV_L1 * (1 + s_`i') 
	gen CF_MW_`i' = BV_L1 * (1 + RonE_`i') - BV_`i' // Clean-surplus dividend: includes repurchases, but treats issuance as negative cash-flow
	replace PV_CF_MW = PV_CF_MW + CF_MW_`i' * (1 + `discount_rate')^(-`i')
	replace t_PV_CF_MW = PV_CF_MW + `i' * CF_MW_`i' * (1 + `discount_rate')^(-`i')
	replace s_L1 = s_`i'
	replace RonE_L1 = RonE_`i'
	replace BV_L1 = BV_`i'
}

gen PV_CF_CS = 0
gen t_PV_CF_CS = 0
* Realized CS-based duration

forvalues i = 1/`horizon' {
	gen CF_CS_`i' = BE_L1_`i' * (1 + ROE_`i') - BE_`i' // Clean-surplus dividend: includes repurchases, but treats issuance as negative cash-flow
	replace PV_CF_CS = PV_CF_CS + CF_CS_`i' * (1 + `discount_rate')^(-`i')
	replace t_PV_CF_CS = PV_CF_CS + `i' * CF_CS_`i' * (1 + `discount_rate')^(-`i')
}


gen w_MW = PV_CF_MW / ME_0
gen Dur_MW = w_MW * t_PV_CF_MW / PV_CF_MW + (1-w_MW) * (`horizon' + (1 + `discount_rate')/`discount_rate')
gen w_CS = PV_CF_CS / ME_0
gen Dur_CS = w_CS * t_PV_CF_CS / PV_CF_CS + (1-w_CS) * (`horizon' + (1 + `discount_rate')/`discount_rate')

bys binq_ME_L1 binq_bm_L1 (year): gen Dur_MW_L1 = Dur_MW[_n-1] if year == year[_n-1]+1
bys binq_ME_L1 binq_bm_L1 (year): gen Dur_CS_L1 = Dur_CS[_n-1] if year == year[_n-1]+1

summarize Dur_MW Dur_CS, detail


********************************************************************************
* 05.05  Generate long-run variables
********************************************************************************	
local rho = .96
foreach var in r roe iva pva ag {
	* Initialize the long-run sum
	qui gen LR_`var' = `var'_0
	forvalues i = 1/19 {
		* Add to the long-run sum
		qui replace LR_`var' = LR_`var' + `rho'^`i' * `var'_`i'
		* Save the long-run sum if a muliple of 5
		local j = `i' + 1
		qui gen LR_`var'`j' = LR_`var' 
	}
}
 
* The remainder: long-run mb
forvalues i = 1(1)19 {
	local j = `i'+1 
	qui gen LR_mb`j' = `rho'^`j' * mb_`i'  
	qui gen LR_Dmb`j' = `rho'^`j'*mb_`i' - mb_L1_0
}


foreach var in r roe iva mb ag pva{
	gen LR_`var'1 = `var'_0
}
replace LR_mb1 = `rho' * LR_mb1
gen LR_Dmb1 = `rho' * mb_0 - mb_L1_0

* Identifier for the combined portfolio bins
gen bin = 10 * binq_ME_L1 + binq_bm_L1 

* Save data
save "Data/Created/NYSE_pfs", replace

********************************************************************************
* 05.06  Cross-sectional mb decomposition 
********************************************************************************

u "Data/Created/NYSE_pfs", clear

* Sample period (first mb_L1_0 available in 1965)
keep if year >= `start' & year <= `end' 
replace LR_mb1 = .96 * LR_mb1

foreach t in 1 5 10 15 20 {
	label var mb_L1_0 " `t' "
	// Change the sign of the effect
	replace LR_r`t' = -LR_r`t'
		foreach var in r roe iva mb {
			eststo: qui reghdfe LR_`var'`t' mb_L1_0 if year <= `end'-`t'+1, absorb(year) cluster(year bin)
		}
		// Save a table for each t then append them togetherk
		esttab, r2		
		esttab using "Results/Tables/05_06_`t'_yr`start'.tex", replace ///
			se(3) nostar r2 b(4) substitute(\_ _) ///
			width(\textwidth) mlabels(none) nonumbers nonotes label ///
			mgroups(" $ -\sum_j=1}^{`t'} r_{t+j} $ " " $ \sum_{j=1}^{`t'} roe_{t+j} $ " ///
				" $ \sum_{j=1}^{`t'} iva_{t+j} $ " " $ \rho^j mb_{t+`t'} $ ", pattern(1 1 1 1)  ///
			prefix(\multicolumn{@span}{c}{) suffix(})   ///
			span erepeat(\cmidrule(lr){@span}))  
		eststo clear
	// Switch the sign back
	replace LR_r`t' = -LR_r`t' 
}

label var mb_L1_0 "$ mb_t $"
	
preserve

// Merge all tables into one 
foreach t in 5 10 15 {
	infix str l 1-500 using "Results/Tables/05_06_`t'_yr`start'.tex", clear
	if `t' < 15 {
		keep if [_n]>=5 & [_n]<=7
	}
	else if `t'==15 {
		keep if inlist([_n],5,6,10,14)
	}
	tempfile tmp`t'
	save "`tmp`t''"
} 
	infix str l 1-500 using "Results/Tables/05_06_1_yr`start'.tex", clear
		keep if [_n]<=7
		append using "`tmp5'"
		append using "`tmp10'"
		append using "`tmp15'"
		replace l = "\hline" if [_n]==2
		replace l = "$ J$ &\multicolumn{1}{c}{ $ -\sum_{j=1}^{J}\rho^{j-1} r_{t+j} $ }&\multicolumn{1}{c}{ $ \sum_{j=1}^{J}\rho^{j-1} roe_{t+j} $ }&\multicolumn{1}{c}{ $ \sum_{j=1}^{J}\rho^{j-1} iva_{t+j} $ }&\multicolumn{1}{c}{ $ \rho^{J} mb_{t+J} $ }\\ " if [_n]==3
		// Replace
		outfile using "Results/Tables/05_06_yr`start'.tex", noquote replace 
	
restore


********************************************************************************
* 05.07  Time-series change in the decomposition (10 year interval) 
********************************************************************************

* Cross-sectional regression

local t = 15 			/* Discounted-sum */
local window = 10		/* Window-regression */
qui sum year, d
local begy = `r(min)' + `window'
local endy = `r(max)' - `t'

* Matrix to store coefficients

foreach var of varlist LR_r`t' LR_mb`t' LR_roe`t' LR_iva`t'{
	mat B`var' = J(1,1,.)
	mat SE`var' = J(1,1,.)
}
	
* Replace the sign of return

replace LR_r`t' = -LR_r`t'

forval i = `begy'/`endy'{
	
	foreach var of varlist LR_r`t' LR_mb`t' LR_roe`t' LR_iva`t'{
		
		* Regressions

		qui reghdfe `var' mb_L1_0 if year < `i' & year >= `i' - `window', ///
			nocon resid a(year) cl(year bin)
		
		* Store coefficients and standard errors
		mat B`var' = (B`var' \ e(b))
		matrix V = vecdiag(e(V))
		matmap V V, map(sqrt(@))
		matrix SE`var' = (SE`var' \ V)	
		
		}
	}
	
* Revert back the sign of return

replace LR_r`t' = -LR_r`t'	

* Restore coefficients and standard errors 

clear
local list "LR_r`t' LR_mb`t' LR_roe`t' LR_iva`t'"
foreach var in `list'{
	svmat B`var'
	svmat SE`var'
}

drop if _n == 1
 
* Set year

g year = `begy'
replace year = year[_n-1] + 1 if _n > 1
tsset year

* Set confidence interval

foreach var in `list' {
	g up_`var' = B`var'1 + 1.96*SE`var'1
	g down_`var' = B`var'1 - 1.96*SE`var'1
}
	
* Label variables

label var BLR_r`t'1  		"LR r"
label var BLR_mb`t'1  		"LR mb"
label var BLR_roe`t'1  		"LR roe"
label var BLR_iva`t'1  		"LR iva"		

* Plot

qui sum year, d
twoway (tsline BLR_roe`t'1 if tin(`r(min)',`r(max)'), lwidth(thick) lcolor(red)) ///
	   (tsline BLR_iva`t'1 if tin(`r(min)',`r(max)'), lpattern(dash_dot) lwidth(thick) lcolor(blue)) ///
		(tsline BLR_r`t'1 if tin(`r(min)',`r(max)'), lpattern(dash) lcolor(gs8)) ///
		(tsline BLR_mb`t'1 if tin(`r(min)',`r(max)'), lpattern(dot) lwidth(thick) lcolor(gs8)), ///
	   ttitle("") graphregion(color(white) lcolor(white)) ///
	   bgcolor(white) ///
	   tlabel(1975 1980 1985 1990 1995 2000 2005) 

gr export "Results/Figures/05_07.pdf", replace as(pdf) 


********************************************************************************
* 05.08 Predicting long-run returns and their cash-flow drivers: with HH std errs (use MATLAB to generate std errs)
********************************************************************************
* Bivariate predictions of long run variables

u "Data/Created/NYSE_pfs", clear

* Sample period (first mb_L1_0 available in 1965)
keep if year >= `start' & year <= `end' 
gen Dur_MW_L1_0 = Dur_MW_L1/10

* Cross-sectionally demean all variables 
local varlist LR_r1 LR_r5 LR_r15 LR_roe1 LR_roe5 LR_roe15 LR_iva1 LR_iva5 LR_iva15 ///
	LR_mb1 LR_mb5 LR_mb15 mb_L1_0 r_L1_0 roe_L1_0 iva_L1_0 ag_L1_0 LR_ag15 Dur_MW_L1_0

foreach var in r roe iva mb ag {
		egen aux = mean(`var'_L1_0), by(year)
		replace `var'_L1_0 = `var'_L1_0 - aux
		drop aux
		forval i=1/20 {
			local j=`i'-1
			egen aux = mean(LR_`var'`i'), by(year)
			replace LR_`var'`i'= LR_`var'`i' - aux
			drop aux
			egen aux = mean(`var'_`j'), by(year)
			replace `var'_`j'= `var'_`j' - aux
			drop aux
		} 
}

foreach var in Dur_MW {
		egen aux = mean(`var'_L1_0), by(year)
		replace `var'_L1_0 = `var'_L1_0 - aux
		drop aux
}
	
* Label regressors		

label var r_L1_0 			" $ r_t $ "
label var mb_L1_0 			" $ mb_t $ "
label var roe_L1_0 			" $ roe_t $ "
label var iva_L1_0 			" $ iva_t $ "
label var ag_L1_0 			" $ ag_t $ "
label var Dur_MW_L1_0 		" $ Dur_t $ "

/* Baseline */
		
* Matrices to store the results
mat B = J(1,6,.)
mat SE = J(1,6,.)
mat Nobs = J(1,1,.)		
		
* Reported results: r, roe, or iva on the LHS
local longrun_var "LR_r LR_roe LR_iva"
foreach lhs in `longrun_var'{
	mat Baux = J(1,1,.)
	mat SEaux = J(1,1,.)
	foreach rhs in r roe iva ag Dur_MW {
		reghdfe `lhs'1 mb_L1_0 `rhs'_L1_0 if year <= `end' - 14, ///
		noabsorb cl(bin year) nocon
		mat Baux = (Baux , _b[`rhs'_L1_0] )
		matrix V = vecdiag(e(V))
		matmap V V, map(sqrt(@))
		matrix SEaux = (SEaux, V[1,2])
	}
	matrix B = (B \ Baux)
	matrix SE = (SE \ SEaux)
	matrix Nobs = (Nobs \ e(N))

	mat Baux = J(1,1,.)
	mat SEaux = J(1,1,.)
	foreach rhs in r roe iva ag Dur_MW {
		reghdfe `lhs'5 mb_L1_0 `rhs'_L1_0 if year <= `end' - 14, ///
		noabsorb cl(bin year) nocon
		mat Baux = (Baux , _b[`rhs'_L1_0] )
		matrix V = vecdiag(e(V))
		matmap V V, map(sqrt(@))
		matrix SEaux = (SEaux, V[1,2])
	}
	matrix B = (B \ Baux)
	matrix SE = (SE \ SEaux)
	matrix Nobs = (Nobs \ e(N))

	mat Baux = J(1,1,.)
	mat SEaux = J(1,1,.)
	foreach rhs in r roe iva ag Dur_MW {
		reghdfe `lhs'15 mb_L1_0 `rhs'_L1_0 if year <= `end' - 14, ///
		noabsorb cl(bin year) nocon
		mat Baux = (Baux , _b[`rhs'_L1_0] )
		matrix V = vecdiag(e(V))
		matmap V V, map(sqrt(@))
		matrix SEaux = (SEaux, V[1,2])
	}
	matrix B = (B \ Baux)
	matrix SE = (SE \ SEaux)
	matrix Nobs = (Nobs \ e(N))
	}
	
tempfile long_horizon
save `long_horizon'


	
* Remove first row / column
mat B = B[2...,2...]
mat Nobs = Nobs[2...,1...]
mat SE = SE[2...,2...]	

* Replace SE with HH SE for long-horions (HHSE from MATLAB)
import delimited "Results/Tables/05_08_HHerrorsBi.csv", clear
mkmat v1 v2 v3 v4 v5 v6

mat hhse2 = v1'
mat hhse3 = v4'
mat hhse4 = SE[4,1...]
mat hhse5 = v2'
mat hhse6 = v5'
mat hhse7 = SE[7,1...]
mat hhse8 = v3'
mat hhse9 = v6'

mat HHSE = SE[1,1...]
mat HHSE = (HHSE \ hhse2)
mat HHSE = (HHSE \ hhse3)
mat HHSE = (HHSE \ hhse4)
mat HHSE = (HHSE \ hhse5)
mat HHSE = (HHSE \ hhse6)
mat HHSE = (HHSE \ hhse7)
mat HHSE = (HHSE \ hhse8)
mat HHSE = (HHSE \ hhse9)

* Export results
clear
svmat B	
svmat HHSE
svmat Nobs

xpose, clear

* Tostring
foreach var of varlist v*{
	qui tostring `var', force replace format("%5.3f")
	}	
	
* Number of obs
foreach var of varlist v*{
	replace `var' = substr(`var', 1, 4) if _n == _N
	}	
	
	
* Put t-stats below coefficients
g order = _n
replace order = order/100 + (_n-5) if _n > 5 
sort order

* Put parantheses around SEs
foreach var of varlist v*{
	replace `var' = "(" + `var' + ")" if _n == 2 | _n == 4 | _n == 6 | _n == 8 | _n == 10
	}

* Add variables names
g var = ""
order var
replace var = " $ r_t $ " if _n == 1
replace var = " $ roe_t $ " if _n == 3
replace var = " $ iva_t $ " if _n == 5
replace var = " $ ag_t $ " if _n == 7
replace var = " $ Dur_t $ " if _n == 9
replace var = " Observations " if _n == 11

* Add linespace after SEs
replace var = " \addlinespace " + var if _n == 3 | _n == 5 | _n == 7 | _n == 9 | _n == 11

* Add \hline before observations
replace var = "\hline \addlinespace Observations" if [_n]==11

* Drop return for the table
drop if _n == 1 | _n == 2

* Change the variables' order
drop order
g order = _n
replace order = 0.1 if _n == 5
replace order = 0.2 if _n == 6
sort order
replace var = " $ ag_t $ " if _n == 1
drop order

* Export it to Latex
listtab * using "Results/Tables/05_08_bi.tex", ///
	rstyle(tabular) replace ///
	head("\begin{tabular*}{\textwidth}{@{\hskip\tabcolsep\extracolsep\fill}l*{9}{c}}" ///
	"\hline" ///
	"&\multicolumn{3}{c}{ $ \sum_{j=1}^{J}\rho^{J-1} r_{t+j} $ }&\multicolumn{3}{c}{ $ \sum_{j=1}^{J}\rho^{J-1} roe_{t+j} $ }&\multicolumn{3}{c}{ $ \sum_{j=1}^{J}\rho^{J-1} iva_{t+j} $ } \\" ///
	"\cmidrule(lr){2-4}\cmidrule(lr){5-7}\cmidrule(lr){8-10}  $ J$ & 1 & 5 & 15 & 1 & 5 & 15 & 1 & 5 & 15  \\ " ///
	"\hline") ///
	foot("\hline" "\end{tabular*}")  
	
	
********************************************************************************
* 05.09 Impulse Response (Use MATLAB to create the figures)
******************************************************************************** 

* Sample period
local start = 1965
local end = 2022

u "Data/created/NYSE_pfs", clear

* Sample period (first mb_L1_0 available in 1965)
keep if year >= `start' & year <= `end' 
gen Dur_MW_L1_0 = Dur_MW_L1/10

* Cross-sectionally demean all variables 
local varlist LR_r1 LR_r5 LR_r15 LR_roe1 LR_roe5 LR_roe15 LR_iva1 LR_iva5 LR_iva15 ///
	LR_mb1 LR_mb5 LR_mb15 mb_L1_0 r_L1_0 roe_L1_0 iva_L1_0 ag_L1_0 LR_ag15 Dur_MW_L1_0 GHpf

foreach var in r roe iva mb ag {
		egen aux = mean(`var'_L1_0), by(year)
		replace `var'_L1_0 = `var'_L1_0 - aux
		drop aux
		forval i=1/20 {
			local j=`i'-1
			egen aux = mean(LR_`var'`i'), by(year)
			replace LR_`var'`i'= LR_`var'`i' - aux
			drop aux
			egen aux = mean(`var'_`j'), by(year)
			replace `var'_`j'= `var'_`j' - aux
			drop aux
		} 
}

	forval i=1/20 {
		local j=`i'-1
		egen aux = mean(LR_Dmb`i'), by(year)
		replace LR_Dmb`i'= LR_Dmb`i' - aux
		drop aux
		gen Dmb_`j' = mb_`j'
		} 

		foreach var in Dur_MW {
		egen aux = mean(`var'_L1_0), by(year)
		replace `var'_L1_0 = `var'_L1_0 - aux
		drop aux
}

* Regressions (first mb_0, r_0, etc. available in 1964) (with year fixed effects) 

keep if year >= `start' & year <= `end'
tempfile data2
save "`data2'"
 
* Bivariate
use "`data2'", clear
mat BiCum = J(1,6,.)
mat Bi = J(1,6,.)
foreach lhs in r roe iva mb{
	forval i = 1/15 {
		local j = `i'-1
			mat BiCumX = J(1,1,.)
			mat BiX = J(1,1,.)
			foreach rhs in r roe iva ag Dur_MW {
				reghdfe LR_`lhs'`i' mb_L1_0 `rhs'_L1_0 if year <= `end' - 14, noabsorb cl(bin year) nocon
				mat BiCumX = (BiCumX, _b[`rhs'_L1_0])
				reghdfe `lhs'_`j' mb_L1_0 `rhs'_L1_0 if year <= `end' - 14, noabsorb cl(bin year) nocon
				mat BiX = (BiX, _b[`rhs'_L1_0])
		}
		mat BiCum = (BiCum \ BiCumX)
		mat Bi = (Bi \ BiX)
	}
}

clear
svmat BiCum

* Save the transition matrix in csv (cumulative)
export excel "Results/Figures/05_09_IRF_pfCUM_Bi.xlsx", replace
clear 
svmat Bi

* Save the transition matrix in csv
export excel "Results/Figures/05_09_IRF_pf_Bi.xlsx", replace

